N <- 750
sims <- 1000

hold <- matrix(NA, nrow = sims, ncol = 3)
set.seed(02140)
for (i in 1:sims) {
  u <- rnorm(N, 0, 2)
  x <- 1.1 * u + runif(N)
  v <- rbinom(N, size = 1, prob = boot::inv.logit(-4 + 3 * x - 1.1 * u))
  d <- rbinom(N, size = 1, prob = boot::inv.logit(-0.5 + 0.5 * x + 0.25 * v))
  
  
  y <- 0.4 * d - 0.4 * d * v - 0.2  * v + -0.2 * x + 0.4 * x * v + rnorm(N, 0, 0.2)

  hold[i,1] <- coef(lm(y ~ d * v + x))["d:v"]
  hold[i,2] <- coef(lm(y ~ d + x, subset = v == 1))["d"] - coef(lm(y ~ d + x, subset = v == 0))["d"]
  hold[i,3] <- coef(lm(v*x ~ d:v + d + v + x))["d:v"]
}

colMeans(hold)
apply(hold, 2, sd)
  cols <- ifelse(d == 1, "indianred", "dodgerblue")

mod1 <- lm(y ~ d * v + x)
mod1.alt <- lm(y ~ d * I(1-v) + x)
mod2 <- lm(y ~ d + x, subset = v == 1)
mod3 <- lm(y ~ d + x, subset = v == 0)
full.mod <- lm(y ~ d*v +x*v)

single <- c(coef(mod1.alt)["d"], coef(mod1)["d"], coef(mod1)["d:v"])
single.ci <- rbind(confint(mod1.alt, "d"), confint(mod1, "d"), confint(mod1, "d:v"))
splits <- c(coef(mod2)["d"], coef(mod3)["d"], coef(full.mod)["d:v"])
splits.ci <- rbind(confint(mod2, "d"), confint(mod3, "d"), confint(full.mod, "d:v"))


cairo_pdf(file = "figures/example-coefs.pdf",
          width = 6, height = 3, family = "Verdana", pointsize = 9)
par(mfrow = c(1,2))
plot(x = NA, y = NA, xlim = c(-1.5, 1.5), ylim = c(0, 4), yaxt = "n", ylab = "", bty = "n", xlab = "Effect of treatment", main = "Single interaction")
abline(v = 0, col = "grey")
abline(v = -0.4, col = "grey", lty = 2)
text(x = -0.40, y = 0.25, "True\ninteraction", col = "grey", pos = 2)
points(x = single, y = 3:1, pch = 19, cex = 0.85)
segments(x0 = single.ci[,1], x1 = single.ci[,2], y0 = 3:1)
text(x = single.ci[,2], y = 3:1, c("Moderator = 1", "Moderator = 0", "Interaction"), pos = 4, cex = 0.85)
plot(x = NA, y = NA, xlim = c(-1.5, 1.5), ylim = c(0, 4), yaxt = "n", ylab = "", bty = "n", xlab = "Effect of treatment", main = "Split samples on moderator")
abline(v = 0, col = "grey")
abline(v = -0.4, col = "grey", lty = 2)
text(x = -0.40, y = 0.25, "True\ninteraction", col = "grey", pos = 2)
points(x = splits, y = 3:1, pch = 19, cex = 0.85)
segments(x0 = splits.ci[,1], x1 = splits.ci[,2], y0 = 3:1)
text(x = splits.ci[,2], y = 3:1, c("Moderator = 1", "Moderator = 0", "Interaction"), pos = 4, cex = 0.85)
dev.off()
